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ABSTRACT 

Most  scientists  who  occasionally  have  to  solve  numerically  a differential 
equation  now  own  a hand  held  programmable  calculator  which  will  very  often  be 
adequate.  Since  hand  held  calculators  are  slow,  there  is  particular  need  to 
keep  the  number  of  function  evaluations  to  a minimum.  At  first  thought,  this 
would  seem  to  rule  out  use  of  Rimge-Kutta  methods,  but  recent  developments, 
such  as  those  by  Fehlberg  (mostly  unknown  except  to  specialists),  may  make  them 
competitive  after  all.  In  the  eurea  of  predictor-corrector  methods,  some  varia- 
tions make  excessive  use  of  memory  locations  for  a hand  held  calculator.  An 
analysis  of  such  matters  is  made  in  order  to  advise  as  to  good  procedures  to 
follow,  including  alerting  the  solver  to  methods  that  are  seldom  taught  in 
nvmerical  emalysis  courses  (where  the  emphasis  is  on  the  use  of  large  fast 
computers ) . 
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significance  and  Explcination 

Until  about  1950,  when  a scientist  had  to  solve  a differential  eqpaation 
numerically  he  had  to  do  it  by  laborious  use  of  a desk  calculator.  Since  that 
time  large  scale  computers  have  revolutionized  the  situation  because  they  can  be 
programmed  to  perform  fast  repetitive  calculations  very  efficiently.  In  recent 
years  hand-held  calculators  have  come  on  the  scene.  Ihey  are  adequate  for  many 
calculations,  in  particular  for  solving  some  kinds  of  differential  equations. 
However  the  factors  that  are  important  for  heind  calculators  are  different  from 
those  for  large  scale  computers.  Hand  calculators  have  only  a limited  number  of 

memory  locations  and  are  comparatively  slow. 

When  solving  differential  equations  on  a hand  held  calculator,  there  is 

particular  need  to  keep  the  number  of  calculations  of  the  function  appearing  in 
the  differential  equation  to  a minimum,  else  the  calculation  will  take  too  long. 
At  first  thought,  this  would  seem  to  rule  out  use  of  what  are  called  Runge-Kutta 
methods,  but  some  recent  developments  may  make  them  competitive  after  all  for 
certain  types  of  problems.  Ihese  types  are  identified,  and  a discussion,  with 
numerical  examples,  is  given  how  best  to  use  the  Runge-Kutta  methods. 

Some  other  methods  which  call  for  much  fewer  calculations  of  the  function 
require  more  memory  locations  than  are  available  on  many  hand  held  calculators. 
There  still  remain  some  methods  which  are  modest  both  in  the  number  of  calcula- 
tions ^md  in  memory  requirements.  How  best  to  use  these  on  a hand  held 
calculator  is  explained,  with  numerical  examples. 


SOLVING  DIFFERENTIAL  EQUATIONS 
ON  A HAND  HELD  PROGRAMMABLE  CALCULATOR 
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Dedicated  to  Prof.  Dr.  Johannes  Weissinger  on  his  65th  birthday. 


1.  Preliminaries . The  present  discussion  is  limited  to  initial  value  ordinary 
differential  equations.  One  wishes  to  solve  a system  of  equations 
(1-1)  y,  = fi(x,y.  ,y. , ....  y ) 

1112  p 

yl  = f-  (x,y  ,y  , . . . , y ) 

2 2 12  p 


being  given  the  values  of  y, .y_,  ...,  y at  x = ; here  y!  denotes 

1 2 p 0 1 

dy^/dx.  We  will  discuss  only  the  special  case 
(1.2)  y*  = f(x,y)  , 

being  given  the  value  of  y at  x = x^;  here  y'  denotes  dy/dx.  The  dis- 
cussion of  (1.2)  can  easily  be  generalized  to  a system  of  equations.  See  Conte 
and  de  Boor,  1972,  pp.  365-366.  Incidentally,  sets  of  higher  order  equations 
can  be  reduced  to  a system  of  first  order  equations,  such  as  (1.1).  See 
Conte  amd  de  Boor,  1972,  p.  365. 

Ibe  person  who  only  occasionally  has  to  solve  numerically  a differential 
equation  may  never  have  had  a course  in  numerical  amalysis,  or  may  have  forgotten 
much  of  it.  So  it  seems  necessary  to  maike  the  present  paper  reasonably  self- 
contadned,  presupposing  little  previous  experience.  Even  if  the  reader  has  pre- 
vious experience,  it  was  presumably  on  a large  fast  computer.  For  the  hand  held 
calculator,  considerations  are  sufficiently  different  that  one  cainnot  rely  too 
much  on  such  previous  experience,  “nie  discussion  to  follow  is  comprehensive 
enough  to  cover  the  important  differences. 

Also,  there  is  disagreement  between  reputable  texts  on  various  points.  An 
adjudication  between  them,  with  explanations,  is  included,  in  spite  of  the  fact 


that  this  adds  to  the  bullc  of  the  present  paper. 
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The  overall  procedure  for  solution  is  the  classic  one,  namely  to  try  step 


by  step  to  calculate  approximations  ^ corresponding  to 


the  values  ...  , where  x^  < < x^  < . , We  tsike  y^ 


the  value  given  for  y at  x = x^. 


We  denote  generically  by  h . Usually  h will  be  taken  the 


same  for  a considerable  succession  of  steps.  A change  in  h is  occasionally 
called  for,  but  not  uncommonly  this  involves  some  travail,  and  in  the  main  one 
tries  to  avoid  it.  On  the  other  hand,  all  authorities  agree  that  at  each  step 
(or  at  least  frequently),  tests  should  be  made  to  see  if  a change  of  step  length 
is  needed,  shorter  to  keep  errors  under  control,  or  longer  to  avoid  unduly  ex- 
tended calculation  that  would  result  from  taking  h smaller  than  necessary. 


2.  The  Euler  method.  If  one  has  proceeded  to  a good  approximation  y^  for 


the  value  of  y at  x = x^  , a Taylor’s  expansion  will  give 

(2.1) 


where  x < ^ < x We  evaluate  y'  by  (1.2),  and  get 

n n+1  n 2 

(2.2)  • 


Unless  one  already  knows  the  solution  of  (1.2),  one  has  no  way  to  calculate 
the  final  term  on  the  right  of  (2.2).  Certainly,  if  h is  taken  to  be  small 
enough,  the  final  term  will  be  quite  small,  amd  the  approximation 


(2.3) 


'n+1  -"n 


y + h f (x  ,y  ) 


n n 


till  give  a sufficiently  accurate  value  for  after  which  one  gets  y^^.2 


by  a similar  formula,  then  yjj+3*  etc.  A problem  with  which  we  must  cope  is 


being  sure  that  we  have  taken  h small  enough  that  repeated  use  of  (2.3)  in- 
stead of  (2.2)  does  not  lead  to  serious  error  in  the  end. 

Suppose  we  wish  to  solve 
(2.4)  y'  « -y. 
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given  that  y = 1 when  x = 0 . Let  us  try  to  approximate  the  value  of 

y when  x = 6.  As  the  answer  for  (2.4)  is  obviously  y = e , we  can  say 

that  for  0 < C ^ 6 we  have  |y"(C)|  < 1.  Hence  the  error  in  (2.3)  will  be 
2 

less  than  h /2.  If  we  taUce  steps  of  consteint  length  h , we  will  require  6/h 

2 

steps  to  get  from  x = 0 to  x = 6.  With  an  error  less  than  h /2  at  each 
step,  and  6/h  steps,  the  final  error  will  add  up  to  less  than 


From  this,  it  is  tempting  to  conclude  that  the  method  is  first  order;  that  is, 
the  overall  error  is  roughly  proportional  to  h . For  exanple,  if  we  decreeise 
h by  a factor  of  2 , we  would  expect  to  decrease  the  error  at  x = 6 by 
approximately  a factor  of  2 . 

Within  bounds,  the  conclusion  is  correct.  However,  the  argument  given 
above  to  support  it  is  fallacious.  To  see  this,  let  us  look  at  a specific 
example.  Take  h = 0.1.  By  (2.3),  we  will  get 

= 0.9. 

This  is  too  small  by  about 

0.0048374. 

2 

(In  accordance  with  (2.2),  this  error  is  close  to  h /2.)  However,  the  value 
of  y at  X = 6 is 
(2.5)  e"®  a 0.0024888. 

Thus  our  final  euiswer  is  less  than  the  error  on  the  first  step.  To  suggest 
that  we  can  approximate  the  overall  error  at  x = 6 by  adding  up  such  errors  as 
shown  above  for  each  of  the  60  steps  required  to  get  from  0 to  6 is  sinqply 
not  sound. 

In  fact,  for  the  equation  (2.4),  use  of  (2.3)  with  h = 0.1  gives 
^n+l  * • 


Using  this  60  times  would  give  an  estimate  for  y at  x = 6 of 
(0.9)®°  = 0.0017970. 

This  is  in  error  by  less  than  28%.  For  a procedure  in  which  the  error  on  the 
first  step  was  almost  twice  the  total  final  answer,  this  is  not  bad. 

By  (2.3),  we  get 


(2.6) 


^n+l  = 


for  the  equation  (2.4).  "niis  gives 


'n+1 


It- •••>'„ 


..)}y_ 


e"^{l 


h ^ h ,, 

1*12- 


e (1 


h . 

3+12  - 


The  factor 

(2.7) 

tends  to  remain  consteuit  eind  close  to  unity  for  |h  ( small,  since  2is  h 
increases  the  first  factor  increases  and  the  second  decreases.  So  we  replace 
(2.7)  by  unity,  getting 


(2.8) 

So  at  X = 6,  we  get 


^ -h.,  h , 


2 n 6/h 


2 


-3h 


a e“^  e“^^ 

For  h * 0.1,  we  get 
y - e"®{0. 74082}. 

Hence , by  the  adsove  emalysis , we  expect  the  calculated  value  to  t)e  too  low 
by  about  26%:  it  was  actually  too  lew  by  aJsout  28%. 
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The  final  formula, 


shows  that  (for  a reasonable  range  of  h ) the  relative  error  is  indeed  of  order 

h . Actually  for  h = 0.1  and  the  equation  (2.4),  we  see  by  (2.2)  that  we  have 

about  0.5%  relative  error  at  each  step.  Accumulating  such  relative  errors  for 
60  steps  could  give  cin  overall  error  of  30%,  which  is  about  what  we  got. 

So  sometimes  it  is  absolute  error  that  one  should  accumulate  from  step 
to  step,  but  other  times  it  can  be  relative  error.  For  programmers  who  try  to 

write  general  purpose  differential  equation  solvers  for  large  computers  with 

error  control  built  in,  the  question  of  how  to  heindle  error  accumulation  poses 
formidedale  difficulties.  See  Hull,  et  al. , 1972,  pp.  607-608.  When  one  is 
solving  on  a hand  held  calculator,  one  sees  the  progress  of  the  solution,  step 
by  step.  In  regions  where  it  is  best  to  reckon  by  accumulating  relative  errors, 
one  can  do  so.  However,  when  it  would  be  better  to  accumulate  absolute  errors 
(for  example,  if  one  is  going  to  pass  through  a zero  value  of  y),  the  change 
to  accumulating  absolute  errors  is  easily  made.  Such  flexibility  is  hard  to 
arrange  in  a preset  program  for  a large  computer. 

By  (2.9) , if  we  should  wish  to  get  to  x = 6 with  0.1%  accuracy,  we  should 
have  to  take  h about  1/3000.  ITius  we  would  require  18,000  steps.  This  would 
not  be  too  bad  on  a large  fast  computer,  but  on  a hauid  held  calculator  it  could 
require  hours,  especially  if  our  equation  to  be  solved  were  more  complicated 
than  (2.4).  So  we  need  something  better  than  the  Euler  method. 

Actually,  there  were  two  reasons  for  including  the  present  section  on  the 
Euler  method.  One  is  to  demonstrate  the  importemce  of  maintaining  flexibility 
about  whether  one  is  accumulating  relative  errors  or  absolute  errors.  With  a 
hand  held  calculator,  such  flexibility  is  easily  maintained.  With  a preset 
program  for  a large  fast  computer,  such  flexibility  is  almost  impossible. 
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Thus,  in  Hull,  et  al.,  1972,  the  decision  was  made  to  lock  one's  self  into 


accumulating  cibsolute  errors  (see  their  p.  607).  To  get  to  x = 6 with  a 
0.1%  error  by  this  scheme  would  require  far  more  than  18,000  steps.  Cn  the  other 
hcuid,  the  large  computers  are  so  very  fast  that  it  is  still  a practical  scheme. 

The  other  reeison  for  mentioning  the  Euler  method  is  that  it  is  useful  for 
"looking  ahead."  Using  a few  steps  with  large  h takes  very  little  time  and 
gives  at  least  a general  idea  of  what  one  might  expect  to  encounter  for  some 
distance  ahead.  If  one  takes  such  a "look  ahead"  every  so  often  it  can  help  in 
planning  how  best  to  arrange  the  upcoming  part  of  the  integration. 


3.  Runge-Kutta  methods.  The  text  Henrici,  1977,  explains  how  to  do  many  calcula- 
tions on  the  HP-25  programmable  calculator.  On  p.  182,  he  suggests  that  for 
solving  the  differential  equation  (1.2)  one  might  use  the  second  order  Runge- 
Kutta  method 

(3.1) 

Appl ied  to 

(3.2)  y'  = k y , 
this  gives 


i'n+l  * ^n  ^ r ^ ^ 


s (l  + hk 


(hk)" 


} 


■'n+1  ^ 2 

An  analysis  like  that  in  the  previous  section  shows  that  if  we  set  y =1  at 


x = 0,  then  at  x = X we  would  get  approximately 

(3.3) 


- ,.-h^k^X/6 

y = e e ' 


'Rie  first  factor  on  the  right  is  what  y should  be,  and  the  second  factor 
shows  about  how  far  off  we  are  from  the  true  value.  So,  if  we  tcOce  k = -1 
(as  in  (2.4))  and  X = 6,  and  ask  for  0.1%  accuracy,  we  need  to  take  h about 
1/32.  Thus  it  will  take  192  steps  to  get  to  x = 6.  However  (note  (3.1)), 
each  step  requires  TWO  evaluations  of  the  function  f(x,y).  So  we  require  384 
function  evaluations. 
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Actually,  we  require  more  than  that.  As  stated  above,  one  should  make 


frequent  checks  to  see  if  h is  about  the  right  size;  this  is  a point  that  we 
failed  to  come  to  grips  with  in  the  previous  section.  Analogous  to  (2.2),  there 
is  a formula  (see  pp.  192-200  of  Ralston,  1965)  that  gives  the  error  of  (3.1)  as 

hK+hL+...  , 

where  K,  L,  etc.,  are  complicated  functions  of  derivatives  of  y and  peurtial 
derivatives  of  f(x,y).  For  small  h , we  may  say  that  the  error  of  (3.1)  is 

3 i 

about  h K.  If  we  take  two  successive  steps,  from  y to  y _ , and  if  K I 

^ •'n  n+2  ' 

does  not  change  much  from  one  step  to  the  next,  then  we  could  accumulate  a local  j 

error  of  2h^K  (over  and  above  whatever  error  we  had  at  y ) in  getting  to  j 

n 

y^^2*  NO'*  3PPly  (3.1)  with  2h  in  place  of  h to  get  from  y^  to  yjj^.2  ^ 

a single  step.  Assuming  that  K is  not  flucuating  badly,  we  would  make  an 

error  of  eibout  (2h)^K  in  getting  to  yjj^.2  ' four  times  the  error  we  made  in 

getting  to  y in  two  steps.  So  now  we  have  two  approximations  for  y _ , 
n+2  n+2 

one  with  an  error  about  four  times  the  other.  From  this,  we  can  estimate  about 

how  far  off  our  better  approximation  is.  If  it  is  close  enough,  we  take  it  as 

the  value  of  y If  it  is  not  close  enough,  we  have  to  go  back  to  y and 

•'n+2  n 

try  over  again  with  a smaller  value  of  h . See  Hull,  et  al.,  1972,  bottom  of 

p.  616.  Hew  close  is  "close  enough"  is  a sticky  question  for  which  there  seems 

not  to  lae  a very  good  answer.  After  all,  this  involves  only  two  steps,  and  one 

heis  to  worry  edaout  the  accumulation  of  errors  over  many  steps.  And,  as  noted  in 

the  previous  section,  there  is  the  question  if  one  should  be  worried  about 

accumulation  of  absolute  errors  or  of  relative  errors.  But  at  least  one  has 

to  have  some  sort  of  estimate  of  the  step  by  step  errors.  j 

i 

As  pointed  out  above,  if  we  take  two  steps  of  length  h each  to  get  from  , 

i 

y to  y _ , we  probably  have  em  estimate  that  is  in  error  by  edaout  one  fourth 
n n+2 

of  what  we  would  get  if  we  went  from  y to  y in  one  step  of  length  2h. 

n n+2 

If  the  error  were  EXACTLY  one  fourth,  we  could  write  a formula  for  the  true 
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value  of  y . It  is  tempting  to  take  this  formula  as  the  value  for  ' 

It  most  likely  qives  a better  estimate  for  y However,  there  is  no  way  to 

n+z 

say  how  muc.  better;  an  estimate  for  the  step  by  step  error  would  then  not  be 
available.  Also,  one  is  left  with  no  improved  value  for  ^ best 

consider  that  this  formula  merely  provides  an  estimate  of  the  error  after  each 
pair  of  steps. 

For  the  pair  of  steps  of  length  h , we  required  two  function  evaluations 

per  step.  For  the  step  of  length  2h  , we  also  require  two  evaluations,  except 

that  the  evaluation  of  f(x,y)  at  x = x^  has  already  been  made.  So,  for  the 

evaluation  of  y , and  y _ , with  an  estimate  of  error  at  x _ , we  require 
n+1  n+2  n+2 

altogether  five  function  evaluations.  So  for  the  192  steps  to  get  from  x = 0 
to  X = 6,  that  is,  96  pairs  of  steps,  we  require  480  function  evaluations. 

While  this  is  a great  improvement  over  the  18,000  function  evaluations  required 
by  the  Euler  method,  it  could  be  rather  time  consuming  if  f(x,y)  is  at  all 
complex. 

So  we  wish  for  something  better.  We  cannot  manage  anything  better  on  the 
HP-25.  It  has  only  50  program  steps,  and  implementation  of  (3.1)  uses  39  of 
these,  leaving  only  11  program  steps  for  the  calculation  of  f(x,y).  In  many 
cases,  11  program  steps  will  not  be  adequate.  So,  if  we  aure  to  do  much  with 
differential  equations,  we  need  a more  capable  calculator  than  the  HP-25. 


4.  Alternatives  to  Runge-Kutta.  Let  us  suppose  that  we  have  a calculator  at 
least  as  capable  as  the  HP-65.  Most  programmable  calculators  now  on  the  market 
are  appreciably  more  capable  than  the  HP-65,  but  it  sufficed  for  the  calculations 
recorded  in  this  paper.  With  it,  one  can  carry  out  higher  order  Runge-Kutta 
methods  than  (3.1).  Suppose  we  use  the  classical  fourth  order  Runge-Kutta. 

See  (5.6-48)  on  p.  200  of  Ralston,  1965,  or  (6.37)  on  p.  338  of  Conte  and  de  Boor, 
1972.  Incidentally,  if  we  have  two  equations 
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(4.1)  y'  = f(x,y,z) 

(4.2)  z'  = g(x,y,z)  , 

then  the  classical  fourth  order  Runge-Kutta  method  consists  of 


(4.3) 

1 

V . = V + — 

(k. 

+ 2k^ 

+ 2k, 

+ k,  ) 

n+l 

'n  6 

1 

2 

3 

4 

(4.4) 

1 

z _ = z + — 

(i( 

+ 

+ 21 

- ^4^ 

where 

n+l 

n 6 

1 

2 

3 

k,  = h f (x  ,y  ,z  ) 
1 n n n 

£i  = h 


k = h f (x  + - , y + — 
Z X\  2.  Ti  z 


Z tA) 

n 2 


S = ^ ^^"n  ^ r ' ^n  ^ ' "n  ^ i-  ' 

k = h f(x  + - , y + — , z + ^ ) 
3 n 2 n 2 n 2 


— 

h 

g (x 

11 

+ — , y 

+ -r-  , z 

+ — 

3 

^ n 

2 ^n 

2 n 

2 

k^ 

__ 

h 

f (x 

+ h , y 

+ k,  , z + 

4 

n 

n 

3 n 

3 

1. 

s 

h 

g (x 

+ h,  y 

+ k,  , z + 

f,)  . 

4 

n 

n 

3 n 

3 

See  Conte  and 

de 

Boor, 

1972 

» PP- 

365-366 

. Extension  to 

such  as  (1.1) 

is 

obvious 

Incidentally , 

the  error 

in  « 

h K for  a suitable  K , and  the  error  in  (4.4)  correspondingly. 

If  we  wish  to  solve  (2.4)  out  to  x = 6 with  0.1%  accuracy,  we  must  take 
h slightly  less  thcin  3/8.  Taking  h = 3/8  is  probably  close  er  utjh,  which 
would  require  16  steps.  There  are  four  function  evaluations  per  step.  To 
accomplish  the  procedure  suggested  above  for  estimating  the  step  by  step  error 
of  ®fter  a pair  of  steps , we  must  also  make  a single  step  of  length  2h 

from  to  Y„^2*  also  takes  four  evaluations,  but  one  has  been  done 

before.  So  wa  require  11  function  evaluations  for  each  pair  of  steps.  So  we 
need  88  function  evaluations. 
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Can  we  do  better?  In  Enright  and  Hull,  1976,  testing  is  reported  of  four 
sorts  of  methods : 

1.  Rimge-Kutta-Fehlberg  methods. 

2.  Rational  extrapolation  methods. 

3.  Adams  type  predictor-corrector  methods. 

4.  Milne  type  predictor-corrector  methods. 

The  Runge-Kutta-Fehlberg  methods  are  Runge-Kutta  type  methods  improved 
according  to  ideas  in  Fehlberg,  1968,  1969,  and  1970, to  give  easier  ways  of 
estimating  step  by  step  errors.  This  is  purported  to  reduce  somewhat  the 
number  of  function  evaluations.  In  view  of  certain  advantages  of  Runge-Kutta 
methods,  the  Fehlberg  improvements  make  the  Runge-Kutta  methods  fairly  competi- 
tive when  the  function  evaluations  can  be  done  quickly.  However  on  p.  949  of 
Enright  and  Hull , 1976,  some  disadvantages  with  the  Fehlberg  methods  are  reported. 

By  phone  (March,  1978)  T.  E.  HiiLl  of  Toronto  informed  the  present  writer  that 
new  improvements  had  eliminated  certain  of  these  disadvantages,  and  that  IMSL 

j 

has  recently  embodied  these  new  improvements  in  its  Runge-Kutta  software  i 

package.  However,  at  best  these  methods  are  competitive  only  when  the  function 
evaluations  can  be  done  quite  quickly. 

The  rational  extrapolation  methods  derive  from  ideas  of  Gragg,  1965,  and 
were  developed  carefully  in  Bulirsch  and  Stoer,  1966,  In  ESiright  and  Hull , 1976, 
these  methods  are  given  a fairly  good  rating  when  the  function  evaluations  can 
be  done  quickly.  They  seem  complex  to  program  for  a hand  held  calculator,  and 
so  we  will  say  no  more  about  them. 

Eiiright  and  Hull , 1976,  give  their  highest  ratings  to  certain  Adams  type 
predictor-corrector  methods.  However,  their  conclusions  do  not  necessarily  hold 
for  hand  held  calculators  because  of  the  very  limited  number  of  memory  locations 
of  hand  held  calculators.  So  we  will  make  a study  of  Adams  type  and  Milne  type 
predictor-corrector  methods  for  hcind  held  calculators. 

I 

i 
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Incidentally,  one  cannot  merely  take  one  of  the  programs  tested  in  Enright 


and  Hull,  1976,  and  put  it  on  his  hand  held  calculator.  In  order  to  accomodate 
the  wide  diversity  of  vagaries  that  can  arise  in  solutions  of  differential 
equations,  these  programs  have  a large  amount  of  "overhead",  and  are  beyond  the 
capabilities  of  present  hand  held  calculators,  besides  being  so  long  as  to  be 
very  time  consuming  at  the  slow  speed  of  hand  held  calculators.  However,  unless 
one  has  a very  large  system  of  equations  to  solve,  one  Cein  proceed  step  by  step 
on  a hand  held  calculator.  One  "looks  ahead"  periodically  as  an  aid  to  planning 


the  calculation,  one  monitors  the  accumulation  of  errors  as  one  goes,  and  one  is 
alert  for  idiosyncracies  that  might  arise.  If  the  latter  do  arise,  one  can  try 


shifting  methods;  if  worse  comes  to  worst,  one  can  try  a change  of  variables, 
or  more  complex  strategems. 


The  classical  predictor-corrector  method  is  that  of  Milne  (see  (28.1)  and 

(28.2)  on  p.  65  of  Milne,  1953,  or  (5.5-12)  on  p.  182  of  Ralston,  1965): 

(4.5)  y^  = y , + ^ (2y'  - y'  , + 2y'  ) 

n+1  n-3  3 n ■'n-1  ■'n-2 


^ntl  = Vl  ^ I ^^‘^n+l'^n+l’  ^ K * K-1^ ' 

This  uses  two  function  evaluations  per  step.  There  is  the  obvious  one  in 

(4.6),  and  when  one  comes  to  use  (4.5)  for  the  next  step,  one  needs  y'  , , which 

n+1 

calls  for  a second  function  evaluation.  Incidentally,  it  is  required  that 

^n+l-i  ■ ^n-i  = i = 0,1,2, 3. 

To  use  this  method,  we  must  have  y ,,  y ,,  y'  V'  V'  stored, 

n-3  n-i  n-2  n-i  n 

cuid  for  the  next  step  we  will  have  to  have  y^_2  ^n'  ^ storage 

requirements  eunount  to  four  values  of  y and  three  of  y' . In  addition,  one 

has  to  carry  the  current  value  of  x . This  uses  up  eight  memory  locations, 

n 

which  for  all  practical  purposes  are  all  there  ^lre  on  the  HP-65.  This  leaves 
no  memory  locations  to  use  in  evaluating  f(x,y).  So  use  of  the  classical 
Milne  predictor-corrector  would  often  not  be  possible  on  the  HP-65.  For 


1 


calculators  with  more  memory  locations,  the  large  requirement  for  memory 
locations  could  preclude  use  of  the  Milne  predictor-corrector  if  one  wishes  to 
solve  a system  of  even  as  few  as  three  equations.  A more  compelling  reason  not 
to  use  the  Milne  method  is  that  it  is  now  known  to  be  unstable  in  certain 
circumstances  which  arise  not  too  infrequently.  We  shall  produce  predictor- 
corrector  methods  that  have  much  more  modest  memory  requirements  eind  are  stable 
besides.  Meernwhile  the  Milne  predictor-corrector  will  be  used  to  illustrate 
typical  features  of  predictor-corrector  methods. 

As  the  name  would  indicate,  a predictor-corrector  method  embodies  a 
predictor  eind  a corrector.  The  predictor  is  (4.5)  and  the  corrector  is 
(4.7) 

which  is  nothing  more  than  Simpson's  rule  for  integrating  y'  approximately. 


= ^n-l  ^ I ^ ^ ^n-l^' 


The  obvious  disadvantage  of  (4.7)  is  that  one  needs  the  value  of  9®^ 

(c)  (c) 

the  value  of  y , , whereas  by  (1.2)  one  needs  the  value  of  y , to  get  the 
n+1  •'n+1  ’ 

corresponding  value  of  y'  , . Actually,  if  h is  reasonably  small  cind  f(x,y) 

n+1 

is  reasonably  well  behaved,  one  can  get  out  of  this  impasse  by  an  iteration 

scheme;  after  making  a guess  for  y’  , , successively  substitute  the  current 

n+1 

(c ) (c ) 

y' into  (4.7)  to  get  a y , and  substitute  y , into  (1.2)  to  get  a better 
n+i  n+i  n+i 

value  for  ^ the  usual  case  that  will  arise,  this  will  converge  to  a 

(c) 

y^^j^  cind  that  satisfy  both  (4.7)  and  (1.2).  Unfortunately,  it  is  prodigal 

with  function  evaluations. 

So  we  adopt  a compromise.  A predictor  is  given,  in  this  case  (4.5),  which 
produces  a reasonably  close  approximation  for  is  substituted  into 

(1.2)  to  get  a guess  for  • This  is  then  substituted  into  the  corrector; 

the  net  result  is  embodied  in  (4.6).  There  the  process  is  stopped.  We  do  not 
have  BiS  much  accuracy  as  we  could  get  with  a few  more  iterations,  but  we  have 
held  the  total  number  of  function  evaluations  to  two  per  step. 
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Since  we  need  four  values  of  y to  use  the  predictor  (and  at  equally  spaced 
values  of  x) , these  four  values  must  somehow  be  obtained  before  we  can  start  to 
use  this  process.  On  pp.  61-64  of  Milne,  1953,  is  given  a scheme  to  get  four 
starting  values.  It  is  now  generally  agreed  that  the  best  way  to  get  started  is 
to  calculate  y^^,  y^ , and  y^  by  Runge-Kutta. 

One  attractive  feature  of  predictor-corrector  methods  is  the  ease  with  which 
one  gets  an  estimate  of  the  step  by  step  error.  The  error  for  (4.5)  is 


45-  h y (O 


and  that  for  (4.7)  is 


1 .5  (v)  , . 

- ^ h y (O, 


(c) 


provided  that  one  has  values  of  y^^j^  and  that  satisfy  both  of  (4.7)  and 


(1.2),  Then  the  error  in  (4.7)  will  be  about 
(4,8) 


-1.  { _ y<P)} 

29  ^'^n+l  l^n+1^  ' 


Actually,  the  value  of  given  by  (4.6)  is  close  enough  to  what  one  would 

get  by  iterating  with  (4.7)  that  the  error  in  from  (4.6)  is  approximately 

(c ) 

what  one  would  get  by  using  y from  (4,6)  in  place  of  y , in  (4.8). 

n +i  n +1 

Since  one  needs  four  consecutive  values  of  y at  equal  step  sizes  for  the 

Milne  method,  a change  of  step  size  (should  it  be  required)  is  not  easy.  If 

one  has  as  many  as  seven  preceding  consecutive  steps  of  equal  length,  one  can 

double  the  step  size  by  pic)ting  every  other  value  from  the  present  and  preceding 

y's.  However  this  would  require  recovering  the  values  of  y ,,  y and 

n-6  n-4 

y^_^  besides  the  values  that  one  usually  stores.  Or  one  could  carry  on  for 

three  more  steps  before  doubling,  being  careful  to  save  the  )cey  values. 

If  one  wishes  to  halve  the  step  size,  one  cem  use 

(4.9)  y = r|5-{45y  + 72y  + lly  + h (-9y ' + 36y'  + 3y'  )} 

i i^o  n n“i  n n-x  n*2 

n-  j 

'“•‘O'  ' 3 ■ ik 

n- 2 
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(see  p.  208  of  Hamming,  1962,  or  (A57)  and  (A58)  on  p.  451  of  Rosser,  1967). 
Unfortunately,  doubling  or  halving  the  step  size  is  often  not  the  most 


I 


I 

'f 


efficient  change  to  make.  For  a change  of  a different  size,  one  can  use  inter- 
polation formulas,  of  which  (4.9)  and  (4.10)  are  samples,  but  it  might  be 
simpler  just  to  make  a fresh  start,  generating  the  next  three  values  of  y 
by  Run  ge  - Kut  ta . 


5.  Adams  predictor-correctors.  For  the  Adams  method  of  order  r , the  predictor 


(which  is  commonly  called  ein  Adams -Bashforth  formula)  is 

= y + h ^ a.  y'  . + K h^^^  (^) 

n+1  ■‘n  1 ■'n-i  p ' 


(5.1) 


and  the  corrector  (which  is  commonly  called  an  Adams -Moulton  formula)  is 


(5.2) 


(c) 

^n+1 


r-1 


= y +hy  B.y',  . 
n . „ 1 n+l-i 


i=0 


The  h that  appears  is 


h = X - X . , 

n+l-i  n-i 

which  is  required  to  be  the  same  for  i = 0,  ...,  r-1. 


(5.3) 


The  error  term  on  the  right  of  (5.2)  is  based  on  the  assumption  that 

(c) , 


y’ 


n+1 


(c) 


As  indicated  with  (4.7),  values  of  ^n+1  satisfy  both  (5.2) 

euid  (5.3)  can  usually  be  found  by  an  iterative  process.  That  is,  if  one 

* (c ) 

chooses  a that  is  near  the  limiting  » ^nd  then  forms 

y;+l  = ^<’‘n+l'  ^n+l^' 

and  substitutes  this  y'  , for  y'  , on  the  right  of  (5,2),  one  will  get  a 

n+1  n+1 


(c ) — 

value  nearer  to  y , them  y 

■^n+1  •'n+1 


Repetition  of  this  converges  to  y 


(c) 

n+1 


However,  this  is  costly  in  function  evaluations,  and  what  is  mostly  done  is  to 
de  fine 


(5.4) 


(t)  = y„  + h Bq  f{x^^^.  y^ll)  + h B^ 


r-1 


'n+1 


i=l 


(t) 


after  which  one  takes  y , = y . This  holds  the  number  of  function 

n+1  n+1 


4 


i 


-14 


evaluations  to  two  per  step.  (One  may  think  of  the  superscript  "t"  as  standing 
for  "traditional.") 

The  requirements  for  memory  locations  are  one  for  y^,  r for 

(i  = 0,...,r-l),  and  one  for  x . Actually,  on  a calculator  like  the  HP-65, 

n 

one  can  reduce  the  memory  requirements  to  one  fewer  by  judicious  use  of  the  stack. 

Suppose  we  eure  at  x , and  have  y , but  have  not  yet  calculated  y*  . Let  us  have 
^ n n n 

y'  stored  at  memory  location  i , for  i = 1,  ...,  r-1,  with  y at  location 
n-i  n 

r , euid  X at  location  r+1.  We  calculate 
n 


r-1 

1 a.  y'  . , 
. 1 n-i 


i=l 

and  store  it  in  location  1,  having  for  i = r-2 , r-3,  ...,  1 successively 

stored  y'  . in  location  i+1.  Now  calculate 
n-i 

y’  = f(x^,y^)  . 
n n n 

Push  an  extra  copy  of  this  up  in  the  stack  to  hold  momentarily.  Now  add  y^ 

to  what  was  stored  in  location  1,  put  the  extra  copy  of  y'  into  location  1, 

n 

and  proceed  with  the  calculation  of  ^ (5.1).  We  now  have  in  storage  or 

in  the  stack  everything  needed  for  calculation  of  by  (5.4). 

In  the  aJaove,  it  is  assumed  that  the  a's  emd  6's  are  part  of  the  program. 
The  a's  and  6's  are  fairly  simple  numbers,  so  that  this  does  not  seem  to 
overload  the  program,  unless  one  has  a very  large  value  of  r , in  which  case 
one  hopes  that  additional  memory  locations  will  tae  available  to  store  the  a's 
and  6's. 

The  Adams  formulas  of  any  order,  with  an  error  term,  can  be  derived  by  use 
of  the  Newton  backward  difference  formula.  Details  are  given  on  pp.  340-342  and 
350-351  of  Conte  2md  de  Boor,  1972,  where  the  Adams  method  of  order  4 is  derived. 
A rather  diffuse  explanation  is  scattered  through  a nixnber  of  pages  of  Milne, 
1953,  but  on  p.  50  is  a tedsle  from  which  the  coefficients  c£ui  be  calculated  up 
to  order  9.  Henrici,  1962,  gives  coefficients  of  the  predictors  on  p.  194  and 


L5- 


r 


of  the  correctors  on  p.  199,  both  up  to  order  6,  Neither  the  tables  of  Milne 

nor  of  Henrici  give  the  error  terms,  but  these  cam  be  derived  easily  by  taking 
r+1 


in  (5.1)  and  (5,2),  which  will  give  the  values  of  and  K^. 


y = X 

Correctors  up  to  order  9,  with  error  terms,  can  be  got  by  a trivial  modifica- 
tion of  formulas  (A2),  (A4) , (A7),  (A12),  (A20),  (A26) , (A34),  and  (A42)  on 
pp.  446-449  of  Rosser,  1967. 

In  order  to  get  started,  we  need  the  r values  Vg'  ^1'  ' ' ' ' ^r  1 ‘ 

For  the  Adams  method  of  order  two,  we  have  the  predictor  and  corrector 

3 


(5.5) 


(5.6) 


(c ) h . • • » h ,,,  / fc.  V 

^n-Hl  = ^n  ^ r <y;tl  ~ n y • 


We  recognize  (5.6)  as  the  trapezoid  rule.  As  noted  for  the  general  case,  we 
define 


h 

n 2 

u (t) 

and  take  y^^^  = y^^^  . 


,(P). 
"n+l'  ^n+l' 


As  y^  is  given,  we  need  to  obtain  a value  only  for  y^^  to  get  started. 

If  one  wishes  to  avoid  use  of  Runge-Kutta,  the  value  of  y^^  can  be  obtained  by 

iterating  with  (5.6).  The  same  applies  if  one  needs  to  make  a restart  after 

chemging  the  length  of  the  step. 

While  the  values  of  C in  the  error  terms  of  (5,5)  and  (5.6)  will  scarcely 

(c ) 

ever  be  the  same,  if  y'"  is  slowly  varying  then  the  error  of  will  be 

about  one  fifth  that  of  Yj^^^  if  we  have  iterated  on  (5.6)  until  both  (5.6) 

cmd  (5.3)  are  satisfied.  In  practice,  the  value  of  Yjj^^  is  close  enough  to 

this  limiting  value  of  Yj^^l  that  one  can  say  that  the  error  of  Yj^^^  is  about 

(p) 

one  fifth  that  of  • in  other  words  the  step  by  step  error  of  the  new 


V , is  about 
■'n+1 


(5.8) 


.L  ( <t)  _ y(p)  , 

6 '^n+l  ^n+1  ’ • 
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If  it  were  not  so,  then  one  is  probably  using  too  large  a value  of  h . 


One  is  tempted  to  try  for  more  accuracy  by  defining 

Qv  (e)  (t)  1 , (t)  (p)  , 

^n+l  = ^n+l  - 6 <yn+l  ' ^n+l  ’ ' 

(e) 

and  then  teiking  to  ^n+1  ' superscript  "e"  to  denote 

"extrapolated".)  This  does  indeed  seem  to  give  more  accuracy.  If  one  would 


_ L (y(c)  _ (p) 

^n+1  6 '^n+1  ^n+1  ’ ' 

(c ) 

with  the  limiting  value  of  » one  would  have  a third  order  method.  See 

Ralston,  1965,  p.  186.  As  y is  close  to  y , one  has  close  to  a third 

n+1  n+1 

order  method.  So  one  should  have  more  accuracy,  but  one  has  no  way  to  tell  how 
much  more;  there  is  no  way  to  estimate  the  step  by  step  error.  Indeed,  as 
we  shall  shortly  show,  niimerical  exaunples  disclose  a somewhat  erratic  behavior 


^n+l  • 

(e) 

The  doctrine  on  use  of  is  not  clear.  In  Rcilston,  1965,  on  p.  186, 

such  use  is  not  favored.  Indeed,  it  is  there  stated  that  it  affects  the  stability 
properties,  and  this  is  also  affirmed  on  p.  210  in  Hamming,  1962.  However,  in  Hamming, 


1962 , there  is  advocacy  of  more  complicated  schemes  which  make  use  of  an  equivalent 
(e) 

of  ^^6  Adams  methods,  which  are  strongly  stable,  the  present  writer 

(e) 

has  found  no  evidence  that  use  of  Y^^j^  causes  stediility  to  deteriorate, 
though  it  does  indeed  for  some  other  types  of  predictor-corrector  methods. 

V 

However,  the  lack  of  any  way  to  estimate  the  step  by  step  error  if  one  uses 
y^^l  seems  a strong  point  against  it. 

Since  we  have  brought  up  the  question  of  stability,  we  should  note  that 
usage  of  the  term  is  not  uniform.  Many  people,  following  Dahlquist,  1956,  say 
that  st2Q>ility  requires  that  all  roots  of  a certain  difference  equation  should 
be  less  them  unity  in  absolute  value.  In  Hamming,  1962,  on  p.  191,  it  is  pointed 
out  that  according  to  that  definition  no  method  for  solving  y'  » ky  could  be 


steible  for  k > 0 . So  Hanuning  introduces  (on  p.  191)  the  notion  of  relative 
stability,  in  which  roots  can  be  larger  than  xinity  in  absolute  value  provided 
they  are  less  in  absolute  value  them  a certain  selected  root.  In  Definition  5.2 
on  p.  176  of  Ralston,  1965,  we  find  Ralston  tedcing  this  as  the  definition  of 
stability;  what  Hamming  calls  "relative  stability"  is  called  "stability"  by 
Ralston.  We  concur  with  Ralston  in  our  use  of  the  word  "stable."  If  a method 
is  stable  in  this  sense , one  can  safely  carry  out  a numerical  integration  of 
indefinite  extent  by  means  of  it. 

It  turns  out  that  to  achieve  stability  for  predictor-corrector  methods, 
one  must  maintain  certain  bounds  for 
(5.10) 


3y 


The  Scune  is  claimed  to  be  true  for  some  Runge-Kutta  methods.  See  Shampine  and 

Watts,  1977,  p.  270.  Writers  of  texts  on  numerical  euialysis  have  been  very  lax 

cibout  calling  this  to  the  attention  of  their  readers,  and  most  people  prolaably 

harlDor  the  illusion  that  Runge-Kutta  methods  are  stable  under  all  circumstances. 

Fortunately,  the  Runge-Kutta' s of  orders  two,  three,  euid  four  given  in  this  paper 

happen  to  be  relatively  stable  (in  the  sense  of  Hamming)  in  all  circumstances, 

and  hence  give  no  trouble  about  stability.  Far  and  away  most  people  who  use  a 

Runge-Kutta  use  oneof  these  three,  which  is  probedsly  why  almost  no  cases  of 

instability  have  laeen  reported  when  using  Runge-Kutta. 

A discussion  of  why  one  needs  to  set  bounds  for  (5.10)  to  avoid  instability 

is  given  in  Hamming,  1962,  cm  pp.  189-190  and  in  Ralston,  1965,  on  pp.  169-178. 

Certainly,  when  one  is  using  predictor-corrector  methods,  one  should  make 

routine  estimates  of  3f(x,y)/3y  to  check  whether  (5.10)  is  remaining  in 

bounds.  A nisnerical  example  of  what  can  happen  if  one  fails  to  do  this  will 

be  given  in  Section  6.  If  one  has  two  values  y and  y + e i and  e is  fairly 

small,  then  we  get  a good  estimate  by 


(5.11) 


(x,y)  ^ f(x,  y+e)  - f(x,  y) 

3y  c 


-18- 


J 


For  the  calculations  proposed,  we  can  take  x = x , , y + e = y » and 

n+1  n+1 

y = y . Thus  we  can  get  our  estimate  by  (5.11)  without  having  to  perform  any 
n+1 

additional  function  evaluations,  since  f (x  ,,  y ) must  be  calculated  to 

n+1  n+1 

provide  the  value  of  required  in  the  predictor  for  the  next  step. 

The  reason  that  y is  close  to  the  y that  one  would  get  by  iter- 
n+1  n+1 

ating  the  corrector,  (5.6),  is  because  y is  already  fairly  close  to  that 

n+1 

; the  use  of  y in  (5.7)  makes  y even  closer  to  y them  y^^]  ■ 
'n+1  n+1  'n+l  n+1  n+1 

If  we  could  find  a still  closer  approximation  tlum  y^^?  to  use  in  (5.7),  we 

n+1 

could  do  still  better.  A closer  approximation  is 


(5.12) 


(P)  + - y^P^, 

.^0.1  f 


'n+1 


n+1  'n+l 


(t) 


since  this  is  Obviously  we  cannot  use  this  until  we  have  calculated 

^n+1  ' However,  with  h fairly  small,  Yj|^^  “ does  not  vary  hugely  from 

step  to  step.  So  if  we  define 


(5.13) 


(a)  (p)  ^ , (t)  (p), 

'n+1  = ^n+l  ^ 'y-  - y-  ) ' 


then  y should  be  closer  to  them  y^^l  . (We  may  consider  the  "a" 

n+1  n+1  n+1 

as  standing  for  "adjusted.")  This  suggestion  is  made  at  the  bottom  of  p.  201  of 
Hamming,  1962,  but  its  use  is  there  discouraged  although  later  a slight  varia- 
tion of  it  is  strongly  advocated.  (We  shall  discuss  this  Hamming  variation  later. ) 
(a) 


If  we  use  y , as  an  improved  predictor,  we  would  define 
n+1 


(5.14) 


(ta) 

^n+1 


"n  * I 'f'-n.l'  'n'l'  * K'  ' 


, (ta)  , • j * ^t) 

and  then  use  ‘or  Y^^j^  instead  of  • 

Since  y^^f^  is  closer  to  y^^l  than  y^^?  or  y^^!  t would  lae  better 

n+1  n+1  n+1 


'n+1 
(a) 

to  define  y by 
n+1 


'n+1 


/r  let  (P^\ 

<5*15)  Y„^l  " Y„^l  + (y„  - ^n  > • 

Actually,  there  is  trouble  getting  started  with  (5.15).  If  one  has  only  y^ 
and  yj^,  there  is  no  way  to  get  (5.15)  cannot  be  used  to  get  Y^*^  ■ 

Me  will  discuss  this  later.  However,  once  one  has  got  started,  we  will  define 
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by  (5.15)  rather  than  (5.13),  and  use  this  in  (5.14). 

n+1 

, (ta)  . , ^ (c)  ..  (t)  . . 

Though  closer  to  than  not  necessarily 

closer  to  the  true  value  of  y , . Indeed,  numerical  results  which  will  be 

•^n+1 

given  shortly  seem  to  support  the  proposition  that  use  of  rather  than 

y^^^  results  in  somewhat  poorer  results  about  half  of  the  time.  However  it 
results  in  appreciably  better  results  the  other  half  of  the  time.  Also,  its 
behavior  is  quite  a bit  less  erratic  than  for  So,  on  the  whole,  it  seems 

a good  idea. 

To  implement  this,  we  would  have  to  allocate  an  extra  memory  location  to 
store  y^^^^-Y  for  use  in  calculating  Y^j^^  next  step.  If  one  is 


(a) 


pressed  for  memory  space,  one  may  have  to  forego  use  of  Also,  there  is 

difficulty  about  calculating  the  first  instance  of  Yq 

given,  and  y^^  estimated  by  some  means,  one  can  follow  the  suggestion  of 

(a)  (p) 

Hamming,  1962,  on  p.  206,  that  one  would  simply  set  y^  = . That  is,  in 

(5.15),  we  ta)ce  y^^^^  - Y^^^  “ 0 . What  this  cimounts  to  is  that  we  use  yJ^^ 
as  an  approximation  for  . It  may  indeed  be  as  good  cui  approximation  for 

y^  as  we  had  earlier  got  for  y^^  . Once  we  have  done  this,  we  can  continue 
on  by  (5.14)  and  (5.15).  Alternatively,  one  can  generate  an  approximation  for 
y^  by  the  same  means  that  we  got  an  approximation  for  y^^  . Then  we  Ceui  get 
y^^*^  by  (5.13),  which  is  close  to  (5.15),  after  which  one  Ceui  use  (5.15)  for 
subsequent  steps . 

It  is  again  tenpting  to  try  for  more  accuracy  by  defining 


(5.16) 


(ea)  (ta)  1 , (ta)  (p), 
^n+l  = ^n+l  - 6 ‘^n+l  ' ^n+l’  ' 

(ea) 


emd  then  taJcing  fo  be  Yjj^j^  • This  has  the  same  adv^mtages  and  dis- 

(e) 

advantages  as  using  y , , Ixit  appears  to  be  less  erratic. 

n+1 

If  we  are  ta)cing  ^ Y^^*\  then  the  two  function  evaluations 

that  would  be  availeUsle  to  put  on  the  right  side  of  (5.11)  would  be  with 
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For  the  second  order  Adams  method  that  we  are 


y + G = y , and  y = y , . 

^ ^n+1  n+1 

considering,  this  should  work  well.  However,  in  Hamming,  1962,  on  p.  207, 

concern  is  expressed  that,  for  a high  order  Adeuns  method,  and  y 

■'n+l  •'n+1 

might  be  so  close  together  that  e and  the  numerator  on  the  right  of  (5.11) 

could  be  appreciably  altered  by  round  off  error,  so  that  use  of  (5.11)  would 

give  a poor  estimate  of  9f(x,y)/9y.  C»ie  should  not  be  working  close  enough  to 

the  bounds  for  (5.10)  that  a sharp  estimate  for  9f(x,y)/9y  is  required. 

However,  one  should  investigate  the  round  off  properties  of  his  calculator  to 

see  if  there  is  danger  of  a really  poor  answer  from  (5.11).  If  there  is,  one 

will  occasionally  have  to  use  an  extra  function  evaluation  to  get  an  extra 

value  from  which  to  estimate  9f(x,y)/9y  by  (5.11). 

To  give  some  assessment  of  the  merits  of  the  preceding  procedures,  we  have 

calculated  the  relative  errors  that  would  result  at  x = 6.  This  has  been  done 

2 

for  the  three  equations  y'  = y , y'  = -y,  cmd  y'  = -2xy  , all  with  y(0)  = 1 

2 

euid  for  three  cases  h = 0.1,  0.2,  and  0.3.  Needless  to  say,  for  y'  = -2xy  , 

2 -1 

the  solution  is  y = (1+x  ) . The  row  labelled  R-K  in  Table  1 gives  the 


result  of  the  Runge-Kutta  of  order  two  that  we  discussed  eairlier.  The  rows 

labelled  "t",  "e",  etc.  refer  to  the  cases  where  y , y^^w  etc.  are  used 

■'n+1  ■'n+1 

for  rows  labelled  "tH"  emd  "eH"  will  be  explained  shortly. 


Reference  to  Table  1 will  verify  some  comments  which  were  made  before.  The 


behavior  of  is  distinctly  erratic,  the  most  extreme  case  being  for 

y'  = y and  h • 0.3,  where  yj^^|  gives  a poorer  result  than  y^^j*  will 

be  noted  that  results  using  poorer  than  yj^^l  fot  y'  = y and 

better  for  y'  * -y,  but  on  the  whole  less  erratic  theui  for  y^^'i  . 

n+1 

In  both  Hamming,  1962,  and  Ralston,  1965,  there  is  advocacy  of  using  y 


rather  than  y , 


where  we  define 


.P 


6.96 

X 

10 

1.47 

X 

10 

•3.73 

X 

10 

3.31 

X 

10 

1.54 

X 

10 

•3.40 

X 

10 

4.47 

X 

lo' 

1.19 

X 

10 

9.89 

X 

10 

2.90 

X 

10 

3.88 

X 

10 

2.98 

X 

lo' 

■2.64 

X 

lo' 

1.85 

X 

lo' 

1.60 

X 

10 

9.47 

X 

10 

4.61 

X 

10 

5.40 

X 

10 

2.81 

X 

10 

9.  35 

X 

10 

3.39 

X 

10 

having  defined 


(5.18) 


y.n  ^n+i’  ^ ■' 


we  take  ^n+1  ' superscript  "H"  stands  for  "Hamming.") 

This  appears  at  the  appropriate  place  in  Table  1.  Hamming  further  proposes 


taking 

(5.19) 


(eH)  (tH)  1 , (tH)  (p), 

=yntl  -6^yn+l  -y--’ 


and  using  y , for  y , , though  Ralston  frowns  on  this  (see  p.  186  of 
n+1  n+1 

Ralston,  1965).  The  results  of  this  also  are  shown  in  Table  1. 


Hamming's  rationale  seems  to  be  as  follows.  Subtracting  y 


from  both 


sides  of  (5.19)  gives 


(5.20) 


(eH)  _ (P)  ^ 5 (tH)  _ (p) 
^n+1  ^n+1  6 ^n+1  ^n+l'  ' 


So  we  can  write  (5.17)  as 


(5.21) 


(H)  (p)  ^ , (eH)  (p), 

^ntl  = ^ntl  ^ - y-  ) • 


If  there  is  very  little  change  of  ~ frcan  n to  n+1,  then  we  are 

very  nearly  taking  'to  be  ^ expressed  on  p.  206  of  Hctmming, 

1972,  we  "mop  up"  the  error  of  y^+i'  Indeed,  later  on  the  same  page  there  is 

(h) 

the  suggestion  that  one  might  just  take  to  be  so  cut  down 

the  number  of  function  evaluations  to  one  per  step.  Presumably  one  still  has  to 

(0H) 

go  through  the  evaluation  of  ^t  each  step,  to  have  it  available  for  use 

in  (5.21)  at  the  next  step,  Isut  for  y’  , in  the  next  predictor  we  would  use 

n+1 

' , (H), 

^^’‘n+l'  ^n+l^  • 

Actually,  one  cannot  quite  hold  it  down  to  one  evaluation  per  step,  since  one 
will  occasionally  need  two  evaluations  in  a step  for  use  in  (5.11).  Incidentally, 
Hamming  does  not  make  an  analysis  of  the  stability  of  this  proceeding.  It  has 
the  disadvantage  that  there  is  no  way  to  estimate  the  step  by  step  error. 


Actually,  if  (5.19)  is  to  give  a third  order  method,  then  should 

(c)  (h^ 

really  be  • For  this,  we  should  like  in  (5.18)  to  be  as  close  as 


possible  to  y , . But  the  definition  (5.17),  being  equivalent  to  (5.21), 
n+1 

is  liable  to  make  y closer  to  y^^^^  than  to  y , and  the  earlier  pro- 
■'n+1  n+1  n+1  ^ 

cedure  of  using  y^^f^  would  seem  preferable.  However,  according  to  Table  1, 
n+1 

y seems  to  do  a shade  better  than  y ^ . The  very  erratic  value  at 
n+1  n+1 

y'  = y,  h = 0.3  for  is  disquieting. 

(t  h) 

If  one  just  stops  at  y^  , and  takes  this  to  be  » instead  of  going 

on  to  (this  is  what  is  proposed  in  Ralston,  1965,  on  p.  189),  it  seems 

better  to  use  the  simpler  y^^f^  . From  what  numerical  evidence  we  have, 

n+1 

y^^f^  is  better  than  y in  half  of  the  cases,  and  worse  in  half. 

To  get  started,  we  calculated  y^^  and  y^  from  the  known  solutions  of 
the  equations,  cind  then  (as  suggested  earlier)  used  (5.13)  to  calculate  • 

If  we  use  these  methods  to  solve  y’  = k y , we  get  stability  in  the  reuiges 
shown  in  Table  2.  This  means  that  for  stability  the  bounds  shown  in  Table  2 must 
be  satisfied  by  (5.10).  Ihe  fact  that  one  has  stability  for  all  positive  values 
is  surprising.  This  fact  is  not  particularly  useful,  since  for  large  hk  the 
step  by  step  errors  would  be  so  large  as  to  render  the  method  of  little  value. 


TABLE  2 

Stability  reinges  for  y'  = ky  for  Adams  second  order. 


-0.6  ^ hk 
-0.7  _<  hk 
-0.7  ^ hk 


Observe  that  in  Table  2 , use  of  extrapolated  values  does  not  diminish  the 
region  of  stability.  If  anything,  the  reverse  is  true. 

In  Section  6 we  will  indicate  how  the  values  in  Table  2 were  derived. 


e 

-0.8  _<  hk  ' 

ea 

-0.7  < hk 

eH 

-0.7  < hk 

1 
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We  now  turn  to  third  order  methods.  We  will  use  the  Rvaige-Kutta 


199  of  Ralston,  1965,  or  25.5.8 


1964 


The  third  order  Adams  is  given  by 


Analogously  to  the  second  order  Adams,  we  set 


TABLE  3 


1 


Relative  error  at  x = 6 for  third  order  Adams. 
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To  get  started,  we  can  estimate  and  by  some  means  and  then  take 

' this  amounts  to  approximating  y^  by  what  was  done 

to  get  the  values  in  Table  3.  Alternatively  one  can  estimate  y^^,  y^,  and  y^ 
by  some  means,  and  then  take  either 

vi"’  - * <y - v'o’ 


yj'-i 


or 


y(H)  = (P)  ^ 2-  (y<t)  _ y(p))  . 

M M 10  '^3  ^3  ' ' 

subsequently  one  uses  either  (5.27)  or  (5.30), 

2 

Some  numerical  results  are  given  in  Table  3.  The  column  headed  y'  = -2xy 
is  quite  erratic.  This  bears  out  the  remark  at  the  top  of  p.  210  in  Hcunming, 
1962,  that  the  corresponding  solution  is  often  troublesome  to  approximate  by 
polynomials. 

For  the  equation  y'  = ky,  we  get  stability  in  the  ranges  shown  in  Table  4. 

These  bounds  should  be  satisfied  by  (5.10).  For  the  case  h = 0.3  for 
2 

y'  = -2xy  , the  bounds  are  not  satisfied  by  (5.10)  for  a region  near  x = 1. 

TABLE  4 

Stability  ranges  for  y'  = ky  for  Adams  third  order  . 


t 

-0.8  £ hk 

e 

-0.9  _<  hk 

ta 

-0.5  ^ hk 

ea 

-0.5  ^ hk 

tH 

-0.5  < hk 

eH 

-0. 5 < hk 

However,  we  got  through  the  region  of  instability  in  two  or  three  steps,  which 
were  not  enough  for  the  instability  to  build  up  appreciably.  For  most  of  the 
range  of  integration,  (5.10)  was  well  within  the  bounds  of  stability. 

As  with  the  second  order  Adams  method,  use  of  extrapolated  values  does  not 
diminish  the  region  of  stability. 
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To  get  some  feeling  how  erratic  some  of  the  values  are  in  Table  3,  note  that 


R-K,  t,  ta,  and  tH  are  supposed  to  be  third  order  methods.  Thus  the  error  for 
h = 0.2  should  be  8 times  that  for  h = 0.1,  and  the  error  for  h = 0.3  should 
be  27  times  that  for  h = 0.1.  Ibis  works  out  reasonably  well  except  for  t . 

As  pointed  out  on  p.  186  of  Ralston,  1965,  use  of 


(5.33) 


^n+1  ^n+1  ” 10  ^^n+1  ” ^n+1  ^ 


should  give  a fourth  order  method.  As  y^^f\  and  y are  close  to 

'^n+1  n+1  •'n+1 

(c ) 

/ use  of  (5.26),  (5.29),  and  (5.32)  means  that  e,  ea,  and  eH  should  be 
close  to  fourth  order  methods.  So  we  look  for  the  error  for  h = 0.2  to  be 


16  times  that  for  h = 0.1,  and  the  error  for  h = 0. 3 to  be  81  times  that  for 


h = 0.1.  None  of  e,  ea,  or  eH  comes  very  close  to  such  behavior.  Nor  can  one 


count  on  a striking  increase  in  accuracy  from  using  extrapolated  values.  In 
two  cases  eH  is  only  barely  better  them  tH,  and  in  one  case  ea  is  poorer 


than  ta  . Considering  that  there  is  no  way  to  estimate  step  by  step  error 


1 = 0.1 

R-K 

4.62  X 10~® 

-5.44  X lo"^ 

-8.10  X lo"^ 

t 

-7.65  X 10“^ 

2.73  X lo"^ 

9.87  X lo"^ 

ta 

-1.39  X lo“^ 

1.61  X lo"^ 

1.02  X lo"® 

tH 

-1.35  X lo”^ 

1.69  X lo”^ 

1.02  X lo"^ 

e 

v£> 

1 

O 

X 

7.50  X lo"^ 

-5.63  X lo'^ 

ea 

-8.10  X lo"^ 

-2.88  X lo"^ 

-5.41  X lo"^ 

eH 

-4.32  X lo"^ 

-2.17  X lo~^ 

-5.40  X lo"^ 

(5.39) 

(ta) 

^n+l 

- K-: 

(5.40) 

(ea) 

^n+l 

, (ta)  19  , (ta)  (p), 

= ^n+l  - - ^n+l^ 

(5.41) 

(H) 

^n+l  = 

(p)  ^ 2^  (tH)  _ (p) 

^n+1  270  ^n  ^n  ’ 

(5.42) 

(tH) 

^n+l 

" IK  - 

(5.43) 

(eH) 

^n+l 

(tH)  19  (tH)  (p), 

= ^n+l  ■ rw  ^^n+l  - ^n+l’  ' 

The  matter  of  getting  started  can  be  handled  as  in  the  second  eind  third 
order  Mams  methods. 

Some  numerical  results  are  given  in  Table  5.  The  extrapolated  values  are 

distinctly  erratic,  as  is  also  However,  behaves  fairly  well. 

The  excellent  results  from  the  Runge-Kutta  are  worth  remarking. 

For  the  equation  y'  = ky,  we  get  stability  in  the  ranges  shown  in  Table  6. 

2 

These  bounds  should  be  satisfied  by  (5.10).  For  h = 0.3  for  y'  = -2xy  , 
this  was  not  the  case  for  a short  interval,  not  long  enough  to  cause  any  trouble. 

TABLE  6 

Stability  ranges  for  y'  = ky  for  Adams  fourth  order. 


-0.6  £ hk 

1 

ea 

-0.4  _<  hk 

eH 

-0.4  < hk 

Use  of  extrapolated  values  does  not  diminish  the  region  of  stability. 


6.  Milne  type  predictor-correctors.  In  Southard  and  Yowell,  1952,  is  given 

^4 


(6.1) 


■ = - 


'n+1 


n n-l  n n~l  o 


(6.2) 


(c)  . h /c  • n ■ . t h (iv) 

^ntl  = ^n  ^ l2  <^y;+l  ^ - 24  y • 

cue  will  recognize  that  the  corrector  is  the  saune  as  for  the  third  order 

Adams,  to  wit  (5.24).  So  the  formula  for  v would  be  the  same  as  (5.25). 

n+1 

We  would  have 


(6.3) 


(e) 


(t) 


(t)  _ . (P), 


, = V - - (v'  ' - V ^') 

'n+1  ^n+l  5 '^n+l  ^n+l'  ' 


and  the  other  formulas  auialogously. 

This  predictor-corrector  has  some  good  features.  It  is  third  order,  but 
has  modest  memory  requirements,  and  is  easy  to  get  started,  or  get  restarted 
after  cheuiging  the  length  of  the  step.  However,  it  behaves  poorly  as  to 
stability,  as  can  be  seen  from  Table  7 and  Tcible  8.  The  fact  that  all  positive 
values  of  h)c  are  excluded  for  ea  euid  eH  is  startling. 

Because  the  range  of  stability  for  y^^]  is  so  limited,  one  finds  (3.10) 

n+1 

2 

outside  that  range  in  the  case  h = 0.3  emd  y'  * -2xy  for  an  extended 
period,  and  the  solution  really  blows  up.  Already  at  x = 3,  one  gets  a negative 
value  for  y , after  which  the  errors  compound  catastrophically.  Before  reaching 
X = 6,  the  calculator  stops  on  an  overflow,  because  the  numbers  are  too  large 
for  its  capacity. 

If  one  had  been  monitoring  (5.10),  this  instability  could  have  been  avoided 
by  using  a smaller  h through  the  critical  region. 

The  analysis  of  the  steUoillty  for  **'*5  ^n+1  ^ worthy  of 

special  attention  by  anyone  interested  in  stability.  Neither  the  usu2d  Dahlquist 
criterion  of  stability  nor  the  Hamming- Ralston  criterion  for  relative  stability 
is  applicable,  and  a special  definition  had  to  be  contrived.  Never  mind  the 
details.  The  reader  can  trust  that  if  (5.10)  stays  too  long  outside  the  bounds 
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TABLE  7 

Relative  error  at  x = 6 for  Southard  and  Yowell. 


y'  = y 

y’  = -y 

II 

to 

to 

R-K 

2.31  X io~^ 

2.71  X 10“^ 

3.31  X 10“^ 

t 

-1.64  X lo“'^ 

-4.18  X 10“'^ 

-1.50  X 10“^ 

ta 

-2.32  X 10“'* 

-2.48  X 10“^ 

-2.01  X 10“^ 

tH 

-2.17  X lo“'* 

-2.73  X 10“^ 

-2.09  X 10“^ 

e 

1.18  X lo“^ 

-1.52  X 10“^ 

-1.82  X 10“'  1 

1 

ea 

unstable 

9.43  X 10“^ 

5.77  X 10“"^ 

eH 

__  1 

unstable 

5.05  X 10“® 

1 

4.16  X 10“^ 

— ! 

R-K 

1.70  X io“^ 

2.35  X 10“^ 

3.03  X lo“‘* 

t 

-8.99  X lo“'* 

-7.16  X 10“^ 

4.54  X 10“^ 

ta 

-1.63  X 10“^ 

-1.75  X 10“^ 

7.17  X 10“^  , 

tH 

-1.45  X lo“^ 

-2.20  X Lo“^ 

-4 

-1.34  X 10  1 

e 

1.70  X lo“'* 

-2.84  X 10“'* 

4.00  X 10“^  1 

ea 

unstable 

-4 

1.79  X 10 

-5  ■ 

4.03  X 10 

eH 

unstable 

1.09  X 10“'* 

1.05  X 10“^ 

R-K 

5.30  X lo“^ 

8.56  X 10“^ 

1.19  X 10“^ 

t 

-2.09  X io“^ 

-2.79  X lo“*- 

unstable 

ta 

-4.66  X lo“^ 

unstable 

1.30  X 10"^  1 

tH 

1 

o 

X 

o 

0 

1 

unstable 

3.45  X 10“^  • 

e 

7.79  X 10“* 

-1.72  X 10“^ 

1.59  X 10"'^  t 

j 

ea 

unstable 

unstable 

unstable  * 

eH 

unstable 

unstable 

1.04  X lo"*- 

TABLE  8 


Stability  ranges  for  y'  = 

ky  for 

Southeurd  and  Yowell. 

t 

-0.31  ^ hk 

e 

-0.7  ^ hk 

ta 

-0.25  ^ hk 

ea 

-0.25  £ hk  < 0 

tH 

-0.28  ^ hk 

eH 

-0.28  ^ hk  < 0 

1 


given  in  Table  8,  there  will  be  trouble;  see  the  case  h = 0.3  for  y'  = -2xy 
for  example. 

In  Hamming,  1959,  on  p.  47,  concern  is  expressed  about  the  stability  of 
Southard  and  Yowell,  and  it  is  proposed  to  remedy  the  situation  by  using  a 
different  predictor 
,(P)  _ 


4 


(6.4) 


^nU  = ^n  ^ ^n-l 


y + 2h(y'  - y • ) + • 

n-2  n n-l  3 


This  appreciably  increases  the  number  of  memory  locations  needed  and  requires 
more  starting  values  to  get  started  or  restcirted,  thus  cancelling  out  the  two 
good  features  of  Southard  and  Yowell.  It  does  appear  to  improve  the  stability; 
for  Yj^^l  it  is  increased  to  -0.5  ^ hk.  However,  this  is  considerably  less 
than  for  the  Adams  method  of  order  three.  If  we  compare  the  Adams  of  order 
three  with  what  results  from  (6.4)  we  find  that  both  are  of  order  three,  both 
need  two  extra  starting  values  to  get  started  or  restarted,  but  the  Adams  is 
considerably  more  stable  and  requires  fewer  memory  locations. 

In  Hamming,  1959,  it  was  proposed  to  render  the  classical  Milne  predictor- 
corrector  stable  by  using  a different  corrector,  ncunely 

the  old  predictor,  (4.5),  was  retained.  The  results  of  this  are  shown  in  Table 

2 

9 euid  TeUsle  10.  In  Table  9 there  is  no  column  y'  = -2xy  because  the  memory 

requirements  exceeded  the  capacity  of  the  HP-65  with  which  the  calculations  of 

this  paper  were  performed.  Although  there  is  a region  of  stability,  it  is  less 

them  for  the  Adams  method  of  comparable  order,  namely  order  four.  All  in  all,  the 

Heunming  method  is  not  a contender  for  use  with  hand  held  calculators. 

From  Ralston,  1965,  cme  would  get  the  impression  that  Hamming's  method  is 
(^K) 

very  superior.  It  is  that  gets  the  stamp  of  approval  on  p.  189  of 

Ralston,  1965.  Actually,  would  be  simpler,  and  gives  alsout  the  same 

n+i 

results.  We  do  not  believe  that  is  enough  superior  to 

n+1  n+i 


to  be 


TABLE  9 


Relative  error  at  x = 6 for  Haniming  , 


h = 0.3 


y'  = y 

y-  = _y 

R-K 

4.62  X lo"^ 

-5.44  X lo"^ 

t 

-8.63  X lo"^ 

4.55  X lo“^ 

za 

-1.68  X lo”^ 

2.10  X lo”^ 

tH 

-1.61  X lo"^ 

2.24  X lo”^ 

e 

4.24  X lo"® 

7.52  X lo"^ 

ea 

-8.10  X 10“^ 

-2.55  X lo"® 

eH 

-4.32  X lo"^ 

-1.92  X 10~^ 

R-K 

6.77  X lo"^ 

-9.46  X lo”^ 

t 

-3.87  X lo“^ 

2,37  X lo“^ 

ta 

-2.07  X lo"** 

2.82  X lo“^ 

tH 

-1.92  X lo”'* 

3,35  X lo"^ 

e 

1.07  X lo“'* 

3.73  X lo"** 

ea 

-8.49  X lo"^ 

-1.15  X lo"'* 

eH 

1.60  X lo“® 

-9.44  X lo“^ 

3.16  X 10 
1.05  X lo" 
-7.53  X lo" 
-6.75  X 10‘ 
6.63  X lo" 
2.32  X lo' 
8.02  X lo' 


-5.21  X 10 
unstable 
6.78  X 10 
1.08  X 10 

5.90  X 10 

-1.10  X 10 

-9.56  X 10 


TABLE  10 

SteUoility  ranges  for  y'  = ky  for  Hamming 
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worth  the  extra  trouble  of  programming  and  the  longer  running  time  that  it 
entails. 

Despite  Ralston's  endorsement  of  Hamming,  Enright  and  Hull , 1976,  give  it  a 

very  low  rating  on  pp.  954-955.  Interestingly  enough,  though  Hamming  invented 

the  method,  and  devotes  a lot  of  space  discussing  it  in  Hamming,  1962,  he 

finally  (in  pp.  206-210)  seems  to  favor  something  related  to  an  Adams  method. 

p.  46  of  Hamming,  1959,  the  stability  for  y is  claimed  to  be  good 

n+1 

down  to  about  -0.65.  As  this  disagrees  with  the  value  given  in  Table  10,  we 
will  justify  the  value  in  Table  10. 

For  Hamming's  method,  we  have 


(6.7) 


(eH)  ^ (tH)  __9_  (tH)  _ p) 
^n+1  ^n+1  121  ^^n+1  ^n+l'  ‘ 


,(P) 


Subtracting  from  both  sides  of  (6.7)  gives 


(6.8) 

So,  by  (6.6) 


(eH)  _ (p)  ^ 112,  (tH)  _ (p) 

^n+1  ^n+1  121  ' ^n+1  ^n+l'  ‘ 


(6.9)  y<«’  = y^PJ  + (y  - y^P’) 

^n+2  ^n+2  '^n+1  ^n+l'  ‘ 

As  the  entries  that  are  accepted  are  those  with  a superscript  (eH) , we 
shall  simplify  the  notation  by  omitting  this.  Ihen  by  (4.5)  euid  (6.9),  we  see 
that  for  the  equation  y'  = ky,  we  have 

yn+2  = ^2  ^ " ^n  ^ 


+ (y, 


K-2  ^ ^ <2y^  - y^_^  + 2y„.2)» 


n+1 


This  simplifies  to 
(H) 


4hk, 


yn+2  - ^n+l  ^ yn-2  ' ^n-B  ^ “ 3y„  + 3y^_^  - 2y^_2>' 
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We  have,  by  (6.5), 


(6.12)  y 


- I - y„,: 


(6.13)  yj'f  ■ J <%.l  - Vl  * “‘%*1  - * ^^-2  ' ^Vs' 


* It>«  ffiy„.i  - I2y„  * I2y„.i  - 8y__.^)). 


By  (6.7) 


i.x  112  (tH)  9 (p) 

^n+2  121  ^n+2  121  ^n+2  ’ 

By  (4.5)  and  (6.13),  this  gives 

(6.15)  y = (126y  - 14y  , + 9y  , 

n+2  121  n+1  n-1  n-i 


+ hk  (150y^^^  - 54y^  + 24y^_^  + 42y^_2  - 42y^_3) 

+ (hk)^(112y^^^  - 168y^  + 168y^_^  - 112y^_2)). 

The  solution  of  this  difference  ec[uation  is 
5 

(6.16)  y = I C.  d"  , 

n . 1 1 ' 

1=1 

where  the  are  the  roots  of  the  equation 

(6.17)  121p^  - (126  + 150hk  + 112(hk)^)p'*  + (54hk  + 168(hk)^)p^ 

+ (14  - 24hk  - 168(hk)^)p^  - (9  + 42hk  - 112(hkf)p  + 42hk  = 0 
If  we  take  hk  = -0.4,  we  get 

(6.18)  121p^  - 83.92p^  + 5.28p^  - 3.28p^  + 25.72p  - 16.8  = 0 . 


The  polynomial  has  the  factors 


(6.19) 

(6.20) 
(6.21) 


p - 0.67053  10701 

2 

P + 0.92764  94248p  + 0.45383  23357 

2 

p - 0,95067  20739p  + 0.45625  70288. 


the  zeros  of  (6.20)  have  absolute  value 
0.67367  07918 


and  the  zeros  of  (6.21)  have  absolute  value 


0.67546  80072  . 


For  large  n , the  powers  of  these  qucintities  will  predominate,  cuid  the  values 


of  will  jump  about  very  erratically. 

For  hk  = -0.39,  the  factor  corresponding  to  (6.19)  will  produce  the 
dominant  p,  and  (6.16)  will  approximate  , as  it  should. 

Amongst  the  disadvantages  of  the  Hamming  method  for  hand  held  calculators 
is  its  excessive  requirement  for  memory  locations.  To  improve  this,  one  could 
chcinge  the  predictor  to 


(6.22) 


'n+1 


9y  + 9y  , + y + h(6y'  + 6y'  , ) + ^ y^^^  (C)  • 
n n-1  n-2  •'n  n-1  10 


When  used  with  the  Hamming  corrector,  (6.5),  this  still  gives  a fourth  order 
method.  However,  it  requires  two  fewer  memory  locations.  Also,  it  requires 
fewer  starting  values  to  get  a start  or  restart,  being  in  this  respect  also 
superior  to  the  Adams  method  of  order  four. 

This  predictor  is  one  of  a set  which,  in  a footnote  on  p.  171  of  Hamming, 
1962,  is  dismissed  as  not  being  worth  consideration.  The  combination  of  (6.22) 
and  (6.5)  does  turn  out  to  have  very  poor  stability  characteristics.  For 
one  must  have  -0.13  < hk,  and  for  y^^f^  one  must  have  -0.19  < hk  < 0.28. 

At  this  point,  we  laecame  disheartened,  2md  did  not  check  out  the  other  cases. 
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7.  Conclusions.  Despite  heroic  efforts  to  remedy  various  faults  of  the 


Milne  type  predictor-correctors,  they  are  still  decidedly  inferior  to  the 
Adams  type.  One  would  prefer  if  the  Adcims  type  results  were  less  erratic 
than  they  are,  but  if  the  user  is  diligent  to  keep  (5.10)  within  bounds,  it 
appears  that  one  can  use  the  Adams  methods  safely,  and  without  running  into 
excessi'^e  calculation  times.  An  occasional  "look  ahead"  to  pleui  the  progress 
of  the  calculation  is  advisable. 

Runge-Kutta  methods  have  many  advantages  and,  if  it  happens  that  f (x,y) 
can  be  evaluated  quickly,  they  should  be  given  serious  consideration. 

Although  the  three  Runge-Kutta's  given  in  this  paper  are  relatively  stable 
under  all  circumstances,  so  that  one  can  use  them  without  any  concern  about 
stability,  one  must  still  monitor  the  size  of  h carefully,  since  if  one 
allows  the  step  by  step  errors  to  become  excessive  the  final  results  can  be 
almost  without  meaning. 
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